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Abstract. 

Redistribution of air masses due to atmospheric circulation causes loading 
deformation of the Earth's crust which can be as large as 20 mm for the ver- 
tical component and 3 mm for horizontal components. Rigorous computa- 
tion of site displacements caused by pressure loading requires knowledge of 
the surface pressure field over the entire Earth surface. A procedure for com- 
puting 3-D displacements of geodetic sites of interest using a 6-hourly pres- 
sure field from the NCEP numerical weather models and the Ponte and Ray 
[2002] model of atmospheric tides is presented. We investigated possible er- 
ror sources and found that the errors of our pressure loading time series are 
below the 15% level. We validated our model by estimating the admittance 
factors of the pressure loading time series using a dataset of 3.5 million VLBI 
observations from 1980 to 2002. The admittance factors averaged over all sites 
are 0.95 ± 0.02 for the vertical displacement and 1.00 ± 0.07 for the hori- 
zontal displacements. For the first time horizontal displacements caused by 
atmospheric pressure loading have been detected. The closeness of these ad- 
mittance factors to unity allows us to conclude that on average our model 
quantitatively agrees with the observations within the error budget of the 
model. At the same time we found that the model is not accurate for sev- 
eral stations which are near a coast or in mountain regions. We conclude that 
our model is suitable for routine data reduction of space geodesy observa- 
tions. 
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1. Introduction 

At the level of precision of modern space geodetic techniques the Earth's crust is not 
static, but deformable. The Earth's crust deformation can be caused by processes inside 
the Earth, by gravitational forces of external celestial bodies, by changes of the centrifugal 
potential, and by various mass loads. Analysis of geodetic observations made from the 
deformable surface of our planet requires applying a model of these deformations. The 
precision of this model should be comparable with the precision of the measurements, 
otherwise unaccounted site position variations due to crust deformation become a fac- 
tor which limits the accuracy of measurements, and the potential of geodetic techniques 
cannot fully be exploited. 

In this paper we focus on the Earth's crust deformation caused by the load of the 
atmosphere. As was found by E. Torricelli in 1644 [Magie, 1963], the atmospheric pressure 
is not constant, but has variations at the level of 20-50 mbar. Darwin [1882] was the first 
who realized that this can cause deformation at the level of several centimeters, and he 
proposed a simple model for its computation. However, before the advent of space geodesy, 
quasi-random site displacements at the level of centimeters were not directly measurable 
and, therefore, there was no necessity to take them into account. 

Rapid development of space geodetic techniques in the 1980-s made it feasible to try 
to detect atmospheric pressure loading signal from the measurements of site positions. 
Trubytsyn and Makalkin [1976] and later Rabel and Zschau [1985], Rahel and Schuh [1986] 
and van Dam and Wahr [1987] made quantitative assessments of the impact of passing 
cyclones and anticyclones on measurements of site positions assuming that the pressure 
distribution in cyclones or anticyclones can be described by a simple mathematical model. 
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Manabe et al. [1991] tried to find a correlation between predicted atmospheric pressure loa- 
ding and the time series of site position determined from very long baseline interferometry 
(VLBl) observations during 1984-1989 at several stations, but had to acknowledge that 
"the observationally determined vertical displacements is not explainable as being caused 
by the atmospheric loading, since the dispersion of the observed vertical displacements is 
too large" . Three years later van Dam and Herring [1994a] and independently MacMillan 
and Gipson [1994] succeeded in detecting atmospheric pressure loading signal using more 
sophisticated approaches. MacMillan and Gipson [1994] estimated coefficients of hnear 
regression between vertical site displacements and local pressure using VLBI data. They 
found that these coefficients for the majority of sites are in reasonable agreement with 
the coefficients derived by Manabe et al. [1991], and that applying an empirical model 
based on the regression between vertical site position and local pressure improves base- 
line length repeatability, van Dam and Herring [1994a], hereafter referred as VDH, used 
another approach. They analyzed the reduction of variance of the estimates of baseline 
lengths derived from analysis of the same VLBI dataset. They found that the reduction of 
variance is consistent with the hypothesis that only approximately 60% of the computed 
pressure loading contribution is present in the VLBI length determination. Applying the 
same technique to global positioning system (GPS) data allowed van Dam et al. [1994b] to 
conclude that 57% of the pressure loading signal is present in the baseline length residuals. 

Although the presence of atmospheric pressure loading signal was confirmed in observa- 
tions, modeling this signal did not come into practice for routine data reduction. First, a 
rigorous computation of displacements caused by mass loading requires handling a gigan- 
tic volume of information and enormous processor power. Second, there was no certainty 
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whether the model is correct. Results of VDH and van Dam et al. [1994b] mean that ob- 
servations did not confirm quantitatively the atmosphere pressure loading model. Without 
solving this discrepancy, applying the atmosphere pressure loading model for processing 
routine observations is not warranted. 

There are several factors that motivated us to re- visit this topic. First, the National 
Centers for Environmental Prediction (NCEP) Reanalysis project [Kalnay et al, 1996] 
now provides a continuous, uniform dataset of surface pressure on a 2°.5 x 2°.5 grid with a 6 
hour resolution for more than 40 years which was not available a decade ago. Second, rapid 
development of high speed networks and processor power makes it possible to retrieve and 
process voluminous meteorological data assimilation models in almost real time. Third, 
the accuracy of geodetic observations has increased considerably during the last ten years, 
which has improved our ability to detect subtle Earth's crust motions. 

The objective of our study was to develop a procedure of computing displacements 
caused by atmospheric pressure loading which is suitable for routine analysis of geodetic 
observations, and to compare these time series of pressure loading with a dataset of 
all VLBI observations from 1980 to 2002. The purpose of this comparison is to get a 
quantitative measure of the agreement between the model and observations, and infer 
whether the model is correct or wrong. In order to do it, we thoroughly examine the error 
budget and on the basis of these estimates compute the expectation of the deviation of the 
observations from the model. Our goal is to determine whether the observations deviate 
from the expectation or not. If the agreement test deviates from the expectation at a 
statistically significant level, this means that either there is an error in computations or 
there is a fundamental flaw in our understanding of the physics of the phenomena under 
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consideration. Then the model must be rejected at this point. If the outcome of the 
statistical test is within the predicted range based on known deficiencies of the model, 
this means that the procedure for computation of the atmosphere pressure loading can be 
accepted for routine reduction of observations. 

In the second section of the paper we describe our method of computing site displace- 
ments caused by pressure loading and assess the error budget of our calculations. We 
re-analyzed a dataset of 3.5 million VLBI observations and performed several statisti- 
cal tests of agreement between the model of pressure loading and the observations. The 
method of data analysis is described in section 3. The results of the VLBI data analy- 
sis are discussed in section 4. Concluding remarks are given in section 5. An efficient 
procedure for computing a time series of pressure loading is outlined in the Appendix. 

2. Computation of Displacements Using Meteorological Models 

According to Farrell [1972] the vertical displacement at a station of coordinates r in- 
duced by surface pressure variations AP(r t) is equal to: 

Ur{f, t) = JJ AP{f', t) Gr{iP) cos ip'dX'dif' (1) 

The vertical Green's function is 

Gn{i') = ^j:h'nPn{cOSi,) (2) 
90 n=0 

f, a and go are the universal constant of gravitation, the mean Earth's radius and the mean 
surface gravity as defined in PREM [Dziewonski and Andersen, 1981], (p' is the geocentric 
latitude, and A' is the longitude. V is the angular distance between the station with 
coordinates r and the pressure source with coordinates f'. Pn is the Legendre polynomial 
of degree n. 
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The horizontal displacement is computed this way: 

Uh{r, i) = JJ ^(r, r ') AP{f', t) Ge(V') cos ip'dX'dip' (3) 

where q{f, r') is the unit vector originating from the station, tangential to the Earth's 
surface, which lies in the plane determined by the radius- vectors to the station and to the 
pressure source. The tangential Green's function is [Farrell, 1972]: 

r f.iA /a+2g,, aPn(cos^) 

Numerical evaluation of the Green's function requires the computation of load Love num- 
bers h'^ and up to a high spherical harmonic degree (n = 9000 in this study) for a 
spherically symmetric, non-rotating, elastic and isotropic (SNREI) Earth model. The 
method of numerical computation of Green's functions is presented by Farrell [1972]. 

We model the oceanic response to atmospheric pressure forcing as an inverted barometer 
(IB): 

AP„ + AP^-AP„ = (5) 

where APa is the variation of local atmosphere pressure, AP^ is the local variation of the 
ocean bottom pressure due to induced sea level change, and APq is the mean atmosphere 
pressure over the world's oceans: 

J J AP{f",t) cos if' dX'd(p' 

APo = (6) 

/ / cos ip' dX' d(p' 

ocean 

which is applied uniformly at the sea floor [van Dam and Wahr, 1987]. This term is 
introduced in equation 5 in order to enforce conservation of ocean mass. Thus, the total 
ocean bottom pressure, AP^ -|- AP^,, is described by equation 6. It has been shown in 

numerous studies (see, for example, [Tierney et al, 2000]) that this model adequately 
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describes the sea height variations for periods longer than 5-20 days. However, the ocean 
response significantly deviates from the IB hypothesis for shorter periods [Wunsch and 
Stammer, 1997]. 

Since APg is zero over the land and depends only on time over the world's oceans, it is 
convenient to split integrals 1 and 3 into a sum of integrals over the ocean and over the 
continental surface. In our computation we use the land-sea mask from the FES99 [Lefevre 
et ai, 2002] ocean tidal model with a 0°.25 spatial resolution. A practical algorithm for the 
computation of site displacements caused by atmospheric pressure loading is presented in 
the Appendix. 

Since the NCEP Reanalysis numerical weather models have a time resolution of 6 hours, 
the semidiurnal {S2) atmospheric tide induced by solar heating cannot be modeled cor- 
rectly, because its frequency corresponds exactly to the Nyquist frequency. The diurnal 
(Si) atmospheric tide is somewhat distorted as well, because of the presence of the ter- 
diurnal signal, which is folded into the diurnal frequency due to sampling. This problem 
was investigated by van den Dool et al. [1997] in detail, who proposed a temporal inter- 
polation procedure which to some extent allows one to overcome the problem. Based on 
this approach Ponte and Ray [2002] recently developed a gridded global model of the 
and 5'2 atmospheric tides with a spatial resolution of 1°.125 x 1°.125. For reasons discussed 
by these authors, we believe their model better represents the atmospheric tides than the 
diurnal and semidiurnal signal which is present in the NCEP Reanalysis model. Using 
these maps we have computed amplitudes and phases of the loading caused by diurnal 
and semidiurnal atmospheric tides for VLBI and sateUite laser ranging (SLR) sites. 
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For each grid point we have estimated four parameters using least squares (LSQ): mean 
pressure, sine and cosine amphtude of the signal, and cosine amplitude of the 5*2 signal 
in the surface pressure field over the time period from 1980 to 2002. This four-parameter 
model is subtracted at each point of the grid from the NCEP Reanalysis pressure field 
before evaluation of the convolution integral. Thus, our time series has zero mean and no 
signal at and ^2 frequencies. The total loading is the sum of the time series and the 
harmonic model of the and <S'2 loading caused by atmospheric tides. 

2.1. Chciracteristics of the Atmospheric Pressure Loading Displacements 

In figures 1 and 2 we show examples of time series of the modeled displacements for the 
period 2000-2003 and their power spectrum at the Wettzell and Hartrao stations respec- 
tively. These stations are representative of mid-latitude and equatorial inland stations. 

All station displacements show significant narrow-band diurnal and semidiurnal signals. 
The displacements for low-latitude stations are characterized by a strong wide-band annual 
and semi-annual signals and relatively weak signal for periods below 10 days, except strong 
peaks at the Si and 5*2 frequencies. Mid-latitude stations show just the opposite behavior. 
For the mid-latitude regions circulation of low and high pressure structures with periods 
of 5-10 days is typical. These periods are also at the edge of the validity of the IB model 
for describing the oceanic response to atmospheric pressure forcing. 

The rms of vertical and horizontal displacements is presented in columns 4 and 5 of 
table 1. Figure 3 shows an example of the temporal autocorrelation function for the 
vertical displacement for the mid-latitude station Algopark. The autocorrelation rapidly 
drops to the level less than 0.2 for time intervals longer than 2 days. The smoothed spatial 
autocorrelation of the vertical displacement induced by atmospheric pressure loading is 
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presented in figure 4. We see that the correlation for baselines shorter than 1000 km is 
very high, typically greater than 0.9, and only for basehnes longer than 3000 km does it 
drop below the level of 0.2 . 

2.2. Error Budget 

There are four major sources of errors in the computation of site displacements caused 
by atmospheric pressure loading: 1) errors in the Green's functions, 2) errors in the 
land-sea mask, 3) errors in the pressure field, and 4) mismodeling the ocean response to 
atmospheric pressure forcing. 

The Green's functions are computed for a SNREI Earth model adopting PREM elastic 
parameters. Thus, we neglect the effects induced by Earth's anelasticity and ellipticity. 
The differences between our Green's functions and Green's functions for an anelastic Earth 
model (see for example Pagiatakis [1990] or Okubo and Tsuji [2001]) are typically below 
1-2%. The effect of the Earth's ellipticity is of the order of magnitude of the Earth's 
flattening, i.e. 0.3%. 

Since the 2°5 spatial resolution of the NCEP Reanalysis surface pressure field is not 
sufficient to correctly represent the coasthne, we chose a land-sea mask with a 0°. 25 res- 
olution from the FES99 [Lefevre et al., 2002] ocean tidal model. The land-sea mask and 
the station distribution are shown in figure 5. We assume that enclosed and small semi- 
enclosed seas (such as Baltic, Black, Red, Caspian Seas and Persian Gulf) respond to 
atmospheric pressure forcing as a "non-inverted barometer", i.e. atmospheric pressure 
variations are fully transmitted to the sea floor. In order to evaluate the effect of the 
errors in the land-sea mask, we computed the time series of atmospheric pressure loading 
with both the 0°.25 and 0°.50 land-sea masks. The differences between these two estimates 
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are typically about 5%. It is worth mentioning that the difference between the loading 
estimates with the 2°. 5 and 0°. 25 land-sea masks can reach 10% for the vertical component 
and 30% for the horizontal components, even for a station like Wettzell (Germany) which 
is 500 km from the coast. 

Another source of error in our computations are errors in the surface pressure field from 
the NCEP Reanalysis. One way to compute this is to compare directly the difference 
between the model and surface pressure observations. Velicogna et al. [2001] presented 
estimates of the rms differences for two different regions (Arabic Peninsula and United 
States) . On the basis of their estimates we conclude that the errors of the NCEP surface 
pressure field on these selected areas are at the level of 5%. 

Another measure of possible errors in the surface pressure field model is the difference 
between two variants of NCEP numerical weather models: the NCEP Reanalysis and the 
NCEP Operational Final Analysis, although these two models are obviously not indepen- 
dent. The NCEP Operational Final Analysis models is an improved version with respect 
to an earlier model NCEP Reanalysis. The improved horizontal (1° instead of 2°.5) and 
vertical (42 layers instead of 28) resolution allows a better modeling of atmospheric dy- 
namics. We computed station displacements due to atmospheric pressure loading using 
the NCEP Operational dataset with a spatial resolution of 1°.0 for the period from April 
2002 to January 2003. The rms of the differences between the vertical and horizontal 
displacements computed with the NCEP Reanalysis and the NCEP Operational data are 
shown in the 6th and 7th columns of table 1. Correlations between the differences and 
the modeled signal are shown in columns 8, 9. The mean error is about 10% for the 
vertical and horizontal components. It is noticeable that the differences are larger for 
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stations enclosed by mountains, for example, Santiago. This is due to the fact that the 
2°5 or even 1°0 spatial resolution of the NCEP Reanalysis and Operational datasets is 
not sufficient to model the topographic variations in mountainous areas and, therefore, 
the surface pressure variations. Our results are similar to conclusions made by Velicogna 
et al. [2001]. The precision of computation of pressure loading displacements is worse for 
mountainous stations. 

In order to evaluate the errors caused by mismodeling of the ocean response, we 
compared ocean bottom pressure variations, as well as the induced loading effects, 
from two runs of the Coupled Large-scale Ice Ocean (CLIO) general circulation mo- 
del [Goosse and Fichelet, 1999]; the first one is forced by atmospheric pressure, sur- 
face winds and heat fluxes [de Viron et al, 2002], and the other one is forced only 
by surface winds and heat fluxes, i.e. assuming an IB response. The differences be- 
tween these two runs can therefore be interpreted as the departure of the ocean re- 
sponse from the IB assumption. We also validated the bottom pressure changes mod- 
eled by CLIO with measurements from the Global Undersea Pressure (GLOUP) data set 
[http://www.pol.ac.uk/psinslh/gloup/gloup.htinl] and found that the CLIO model 
agrees with the measurements of the bottom pressure at the level of 20%. We show in 
the 10th and 11th columns of table 1 the ratio of the rms of these differences to the rms 
of our atmospheric loading series with the IB model for all stations, as well as a mean 
value of this ratio. Correlations between the differences and modeled signal are shown in 
columns 12 and 13. Thus, the mean vertical and horizontal errors are below 10% and 20% 
respectively. As expected, these values are higher for island stations (f.e., Kokee, Mk- 
vlba) or stations close to the coasts (Richmond, Hobart26, etc.), where the atmospheric 
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loading itself is very small: rms below 1 mm and 0.5 mm for the vertical and horizontal 
components respectively. 

Table 2 summarizes the error budget. Combining all known sources of errors we evaluate 
the total uncertainty of our computation of site displacements due to atmospheric pressure 
loading to be 15%. 

3. Validation of the Model Using VLBI Observations 

We selected VLBI for validation of our time series of atmospheric pressure loading. 
Each of the three main space geodetic techniques, GPS, SLR and VLBI, has its own 
advantages and disadvantages, although in general they are quite competitive. We chose 
VLBI because of the maturity of the VLBI data analysis technique. Complete re-analysis 
of the whole set of VLBI observations takes about a couple of hours on rather a modest 
computer. Therefore, the consistency of reduction models and parameter estimation can 
easily be enforced. These factors make VLBI attractive for investigating tiny effects hke 
atmospheric pressure loading. 

3.1. Observations 

All dual-band Mark-3/Mark-4/K-4 VLBI observations carried out under various 
geodetic and astrometric programs from 1979 to the present are available on-line at 
the International VLBI Service for Geodesy and Astrometry (IVS) Data Center at 
http://ivscc.gsfc.nasa.gov [Vandenberg, 1999]. The VLBI data set has substantial 
spatial and time inhomogeneity. Typically, observations are made in sessions with a dura- 
tion of about 24 hours. Observations were sporadic in the early 80s, but in January 1984 
a regular VLBI campaign for the determination of EOP started first with 5-day intervals. 
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from May 1993 with weekly intervals, and from 1997 twice per week. In addition to the 
observations dedicated to EOP measurements, various other observing campaigns were 
running. On average 150 sessions per year have been observed since 1984. 

In total 144 stations participated in observations, although a majority of them observed 
only during short campaigns. The stations which participated in more than 20 000 ob- 
servations for more than three years were used for analysis. 46 stations satisfied these 
criteria. Four stations which participated in the Key Stone Project (KSP) [Takahashi 
et al., 2000], Kashimll, Koganei, Miura, Tateyama, were excluded since they observed 
mainly in a small local network, as well as two other stations, Crimea, because it had 
bad performance, and Ylow7296, since its sensitivity was too low. Only observations on 
the baselines between the 40 strong selected stations were used, and other observations 
(~6%) were discarded. Sessions with less than 3 strong stations were discarded entirely. 
3073 sessions from April 1980 to December 2002 with more than 3.5 million observations 
remained, and they were used in the analysis. 

The number of participating stations in each individual session varies from 2 to 20, 
although 4-7 is a typical number. No station participated in all sessions, but every 
station participated in sessions with many different networks. All networks have common 
nodes and, therefore, are tied together. Networks vary significantly, but more than 70% 
of them have a size exceeding the Earth's radius. 

3.2. Choice of Parameterization 

The scatter of daily estimates of site positions is greater by a factor of 2-5 than the 
rms of atmospheric pressure loading displacements. Therefore, we cannot directly see the 
signal in the site position time series. Besides, in order to resolve rank deficiency of the 
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problem of simultaneous adjustment of coordinates of all site and EOP, net-rotation and 
net-translation constraints should be applied or an equivalent technique should be used. 
As a result adjustments of site positions become linearly dependent. Since only a small 
number of stations participates in a typical VLBI experiment, the influence of position 
variations of other observing stations of the network on position variations of a station of 
interest is not diluted to a negligible level. It makes the interpretation of daily VLBI site 
position time series uncertain. 

One of the ways to assess validity of the model is to compute two time series of basehne 
lengths: the first with applying the atmospheric pressure loading model and the second 
without applying the model. Baseline lengths are invariant with respect to a rotation 
and translation and, therefore, net-translation and net-rotation constraints do not affect 
them. We introduce the reduction of variance coefficient R: 



where Acr^ is the difference between the mean square of baseline length residuals before 
and after applying the model, and am is the variance of the signal in the model. If 
the model is perfect and the signal under consideration is not correlated with another 
unmodeled effects, the coefficient is 1. If the baseline length series does not contain the 
signal coherent with the model at all, then applying the model increases the variance by 
the amount of the variance of the signal in that model, and the coefficient of reduction of 
variance is 0. We should emphasize importance of the assumption of the lack of correlation 
between the pressure loading signal and noise: we can extract the signal which is below 
the noise level if, and only if, some additional information about the noise is exploited. 
The validity of this assumption is based on the fact that as can bee seen in figures 1 




(7) 
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and 2, the spectrum of the atmosphere pressure loading at frequencies below one day is 
flat, except for a peak at the annual frequency for some stations. Thus, the atmospheric 
pressure loading can be considered to some degree as a stochastic Gaussian process. The 
spectrum of the VLBI residuals is also flat. Therefore, our condition of lack of correlation 
between the atmosphere pressure loading series and residual unmodeled effects is fulfilled 
if an unmodeled contribution to VLBI delay and pressure loading are independent. 

Although this approach gives us a quantitative measure of the adequacy of the model, 
it has some disadvantages. It lets us determine only the reduction of variance coefficients 
for the projection of the difference of the site displacements vectors on the baseline vector 
instead of the coefficients for each site and each component independently. 

Another approach is to represent the atmospheric pressure loading signal as a product 
A ■ ttm, where is the modeled signal, and to estimate directly from the VLBI time 
delays the unknown parameter A which hereafter we call admittance factor between the 
modeled signal and the observables. It can easily be verified that if A is the only estimated 
parameter and the modeled signal is not correlated with an unmodeled contribution to 
the observable, then the expectation of the LSQ estimate of the admittance factor, E{A), 
is 

E{A) = (8) 

where p is the correlation coefficient between modeled and true atmospheric pressure 
loading, and ai is the variance of the true pressure loading. Under the same assumption 
the reduction of variance Aa^ = (^E[A) — l^o"^ and therefore, the coefficient of reduction 
of variance R is equal to E{A). If other parameters are adjusted in addition to A, this 
property is preserved at the level of correlations between A and the estimates of other 



DRAFT 



February 2, 2008, 3:39am 



DRAFT 



PETROV AND BOY: ATMOSPHERIC PRESSURE LOADING X - 17 

parameters. The admittance factor A shows how much of the power of the modeled signal 
is present in the observables. 

3.3. Estimation Model 

We made solutions of two classes: global solutions, in which we estimated the positions 
and velocities of all sites over the entire data set, and baseline solutions, in which we 
estimated site positions for each VLBI experiment independently using ionosphere free 
hnear combination of group delays at S and X bands. Estimated parameters were split into 
two classes: basic parameters, which are usually adjusted in processing VLBI experiments, 
and specific parameters of interest. Basic parameters belong to one of the three groups: 

— global (over the entire data set): positions of 511 primary sources and proper mo- 
tions of 79 sources. 

— local (over each session): pole coordinates and their rates, UTl, UTl rate; positions 
of other sources, azimuthal troposphere path delay gradients for all stations and their 
rates, station-dependent clocks modeled by second order polynomials, baseline-dependent 
clock offsets. 

— segmented (over 0.33-1.0 hours): coefficients of a linear spline which models at- 
mospheric path delay (0.33 hour segment) and clock function (1 hour segment) for each 
station. The estimates of clock function absorb uncalibrated instrumental delays in the 
data acquisition system. 

The rate of change of atmospheric path delay and clock function between two ad- 
jacent segments was constrained to zero with weights reciprocal to 40psec/hour and 
2 • 10~^^ sec/sec respectively in order to stabilize solutions. No-net-translation and no- 
net-rotation constraints were imposed on the adjustments of site positions and velocities 
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as well as no-net-rotation constraints were imposed on adjustments of source positions in 
order to solve the LSQ problem of incomplete rank. 

A pair of baseline solutions was made with only basic parameters estimated and station 
positions treated as local parameters. In reference solution Bl the atmospheric pressure 
loading time series was not applied; in solution B2 the contribution due to the model of 
atmospheric pressure loading was added to the theoretical model of the observables. 

In solutions Gl, G5, G6 and G7 the list of global parameters included site positions, 
site velocities and overall admittance factors for the Up, East and North components of 
the modeled atmospheric pressure loading signal for all sites combined. In solution G2 we 
estimated the set of the admittance factors for the Up, East and North components for 
each station independently, and other parameters were the same as in the Gl solution. 

It should be noted that unlike estimation of positions of all stations, estimation of 
the Up, East and North admittance factors for all stations of a global network does not 
result in a rank deficiency. A partial derivative of time delay r with respect to the 

dT dT 

admittance factor of the i'th station is given by — — ami, where n is a vector of 
station coordinates. As was shown in figure 3, the correlation between the atmosphere 
pressure loading at different stations is less than 0.1 at distances greater than 4000 km. 

3.4. Theoretical Model 

The computation of theoretical time delays, with some exceptions, generally follows the 
procedure outlined in the lERS Conventions [McCarthy, 1996] and described in more detail 
by Sovers et al. [1998]. The GOTOO model [Ray, 1999] of diurnal and semidiurnal ocean 
tides, the NA099 model [Matsumoto et al., 2000] of ocean zonal tides, the equilibrium 
model of the pole tide and the tide with period of 18.6 years were used for the computation 
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of displacements due to ocean loading. The hydrology model of Milly and Shmakin [2002a] 
was used for the computation of displacements due to hydrology loading (D.S. MacMillan 
and J.-P. Boy, Effect of mass loading signals on crustal displacements measured by VLBI, 
manuscript in preparation; D.S MacMillan, personal communication 2003). The empirical 
model of high frequency Earth orientation parameters derived from analysis of the VLBI 
data [Petrov, 2000a] and the IERS96 semi-empirical nutation expansion were used. The 
Niell [1996] mapping functions were used for modeling and estimating the tropospheric 
path delay. The displacement of the reference point of a VLBI station due to thermal 
expansion of the antenna was not modeled. Although Nothnagel et al. [1995] proposed 
a model for an antenna's thermal expansion, attempts to validate this model were not 
successful. We did not include the model of non-tideal ocean loading, because this effect 
is one order of magnitude smaller than the atmosphere pressure loading, and it has not 
been fully investigated. It will be considered in a future paper. 

4. Discussion 

4.1. Analysis of Global Admittance Factors 

Table 3 shows admittance factors determined in solution Gl. The uncertainties of the 
results were derived by propagating the group delay errors. These errors were computed 
on the basis of the signal to noise ratio of the cross correlation function of the recorded 
signal from the receivers and the empirical baseline-dependent reweighting parameters 
which, being added in quadrature to the uncertainties of group delays, made the x^f of 
the postfit residuals close to unity. In addition to that, the formal uncertainties were 
scaled by a factor of 1.5. Numerous tests of splitting the dataset into subsets in time and 
in space showed that the VLBI formal uncertainties are underestimated by the factor of 
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1.5. The origin of this discrepancy is not completely understood, but there are indications 
that it may be caused by fringe phase variations due to instrumental errors [Petrov, 2000b] . 

Above we showed the estimates of the rms of possible errors in our computation of 
atmospheric pressure loading and their correlations with the modeled signal. Assuming 
that the different error sources a) are not correlated with each other, b) are small with 
respect to the signal, we can present the expectation of the admittance factor in this form: 

E{A) ^1 -^kin + J2kl {2r^ - 1) + '^J2Jlk^kJnr, + 0{e) (9) 

where k is the ratio of the rms of errors to the rms of the modeled signal, r is the correlation 
between the error and modeled signal. The summing is done over all considered sources 
of errors. 

Using numerical values for and ki listed in table 2, we evaluate the expectation of 
A: 0.90 . Our estimates of the admittance factors are close to this value. This means 
that the known deficiency of the model is sufficient to explain the small deviation of the 
estimates of the global admittance factors from 1. 

Since the typical spectrum of atmospheric pressure loading shows peaks at semi-diurnal 
(1.46-10"'^ rad/s), diurnal (7.29 -10"^), semi-annual (3.98-10"'' rad/s) and annual fre- 
quencies (1.99 -10"^ rad/s) (figures 1-2), we would hke to see how applying the atmo- 
spheric pressure loading model affects residual harmonic site position variations at these 
frequencies. In order to assess this effect, we made two solutions, G3 and G4, and es- 
timated the sine and cosine amplitudes of position variations of all sites at 32 tidal fre- 
quencies, including these four frequencies with noticeable atmospheric pressure loading 
signal. Atmospheric pressure loading was applied in solution G4, but was not in solution 
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G3. Amplitudes at the frequencies of diurnal, semidiurnal and long period bands where 
no tidal signal is expected were estimated as well in order to calibrate the uncertainties 
of the results. The ratio of the weighted sum of squares of the residual amplitudes over 
all stations at the specific frequency to its mathematical expectation, was used as a 
measure of the power of the residual signal. In the absence of the signal these statistics 
should be less than 1.25 at the 95% confidence level. Therefore, large values of V which 
exceed this limit indicate the presence of the residual signal. This technique is explained 
in more details by Petrov and Ma [2002]. 

Table 4 shows the estimates of V for each of the four frequencies of interest. We 
see that applying the model of atmospheric pressure loading reduces the amplitude of 
the residual signal at semi-diurnal and semi-annual frequencies, but noticeably increases 
the amplitude at the annual frequency. In all cases the power of remaining residual 
signal is still significant. This table shows us that the atmosphere pressure loading is 
not the dominating source of observed residual harmonic site position variations at these 
frequencies. 

The presence of the narrow-band residual signal to some degree violates the assumption 
of lack of correlations which was put at the basis of the agreement test. Two harmonic 
signals may be independent and correlated. In order to investigate the effect of the 
atmospheric pressure loading signal at the annual frequency, the frequency with the largest 
harmonic signal, we passed the pressure loading time series through a narrow-band filter 
with a passing window around the annual frequency: [1.79 • 10^^, 2.19 • 10^^] rad/s. We 
made two solutions, G5 in which we estimated the global admittance factor of the annual 
constituent of the atmospheric pressure loading and solution G6 in which we estimated 
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the global admittance factor with the signal at the annual band removed. The results 
are presented in table 3. The admittance factor in solution G5 increased by 3%. It gives 
us a measure of the distortion of the statistical tests caused by unaccounted correlations 
between the annual narrow-band unmodeled signal and atmospheric pressure loading. 
Similar results were reported by VDH who noted that removal of the annual signal from 
the atmospheric pressure loading increased variance reduction of the baseline length series. 

The fact that modeling atmospheric pressure loading increases the level of the residual 
signal at annual frequency does not necessarily mean that the model is wrong. If two 
anticorrelated signals contribute to the observables, then including the model of only 
one of the signals in the procedure of reduction of observations may result in increasing 
variance, even if the model is perfect. We know that various phenomena may contribute 
to the annual site position variation, such as hydrological signal, non-tidal ocean loading, 
thermal antenna deformation, mismodeling troposphere path delay, etc, and atmospheric 
pressure loading is not the greatest contributor. It should be noted that this means that 
a reduction of the variance coefficient and an admittance factor cannot be considered as 
a valid test of goodness of the model in this specific case. 

We should acknowledge that currently we are unable to test directly whether the annual 
constituent of the atmospheric pressure loading signal is modeled correctly or not. It will 
be possible in the future when a complete model of site position variations will be built 
and the power of the residual signal will become less than the power of annual pressure 
loading signal. Although attempts to model seasonal effects show certain improvement 
[Dong et al., 2002], we are still far from a solution of this problem. 
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At the same time the admittance factor for the vertical displacements is close to unity at 
the level of measurement noise for the wide-band component of the modeled displacements 
due to atmosphere pressure loading after subtraction of the annual component. It provides 
us indirect evidence that we have modeled atmosphere pressure loading correctly at the 
annual frequency as well, since the Green's function and land-sea mask are frequency 
independent, and our estimate of the error budget set the upper limit of possible seasonal 
errors of the atmosphere pressure field. 

4.2. Analysis of Admittance Factors for Individual Stations. 

Although we concluded in the previous section that the average admittance factor is 
very close to unity, it does not necessarily mean that the admittance factors are close to 
unity for each individual station. Table 5 shows the estimates of the admittance factors 
from the G2 solution. The estimates with the formal uncertainties greater than 0.5 are 
omitted in the table. It is not surprising that the admittance factor of the Kokee station 
is far from unity, even negative, since the station is located on an island in the middle of 
the Pacific ocean. It is more surprising that the admittance factors at the Hn-vlba and 
Westford stations are so different, since the distance between these stations is only 54 km 
and, indeed, the time series of the pressure loading are very similar. We suspect some 
problem with data at the station of Hn-vlba. Low admittance at the La-vlba and Pietown 
stations located in mountainous regions may be explained by greater errors in the model 
of surface pressure. We do not have explanation of the negative admittance factor at the 
station of Ov-vlba. Another peculiarity is the anomalously high admittance factors at 
the Seshan25 and Tsukub32 stations. It is known that the response of the Yellow Sea 
to tidal forcing is amplified by 10 times, and a significant non-linear M4 tide is observed 
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[Lefevre et ai, 2000]. We can expect that the response to atmospheric forcing will be also 
substantially non-linear and not consistent with the IB hypothesis. 

4.3. Analysis of Reduction of Variance Coefficients 

In order to compare our results with the early VDH paper, we performed our analysis 
in a manner similar to that used in the analysis of those authors. We computed the time 
series of baseline lengths in the Bl and B2 solutions. A linear model was fit in the series 
with discontinuities at epochs of seismic events for several stations. The weighted root 
mean square of residual baseline lengths was computed for all baselines with more than 
100 sessions for both the Bl and B2 solutions. 69 baselines satisfy this criterion. The 
coefficients of reduction of variance were computed using baseline length variances. 

The histogram of the distribution of the coefficients of reduction of variance is presented 
in figure 6. The weighted mean value R over 69 baselines is 0.97 ± 0.04 . For the com- 
putation of the uncertainty of the mean value we used the variance of R which can easily 
be derived from the results presented in the appendix of VDH: 

^2 I 2 ^2 

var OR = —-^ 10 

where and o"| are the mean squares of the baseline length residuals of solutions B 1 and 
B2 respectively, cr^ is the mean square of the modeled basehne changes with respect to 
the mean value computed over the period of time of VLBI observation, N is the number 
of observing sessions. 

VDH analyzed 22 baselines for the period 1979 to 1992. They reported a much lower 
value of the reduction of variance coefficient, 0.62, which means that 38% of the power 
of the signal in their series of atmospheric pressure loading is not present in the VLBI 
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data. In order to check whether the differences may be due to changes in the quahty 
of VLBI data collected after 1992 we restricted our calculations to exactly the same set 
of observations and baselines used by VDH. We got 1.10 ± 0.10 for the coefficient of 
reduction of variance for this case. It deviates at the 3a level from the value reported by 
VDH. The differences in our analysis technique of VLBI observations and the technique 
of VDH are not significant enough to explain this large discrepancy. So, the differences 
in the reduction of variance coefficients come from the differences between the series of 
the atmospheric pressure loading displacements. First, the predicted baseline loading 
effects of VDH and van Dam et al. [1994b] suffered from the problem that they have the 
wrong sign for the horizontal deformation in the North direction (T. van Dam personal 
communication, 2003). Second, we used the NCEP Reanalysis model, but VDH used 
an older model: the National Meteorological Center operational model with a 12 hour 
resolution. Third, we used a land-sea mask with a resolution of 0°.25 x 0°.25 — ten times 
better then the resolution of the surface pressure model. As shown in section 2.2, this 
difference alone can cause an error in the vertical displacements as large as 10%. Fourth, 
we used a different numerical algorithm for computing the convolution integral. 

4.4. Application of the Pressure Loading Model to Data Reduction 

In the past, several authors recommended finding the linear regression of the vertical 
atmospheric pressure loading and local surface pressure and using this simple model in op- 
erational data analysis. Rabel and Zschau [1985] warned that since "the magnitude of the 
displacements is critically dependent on the spatial extension of the pressure distribution, 
[. . . ] there cannot be any unique regression coefficient between local displacements and 
local air pressure changes which could be used to correct geodynamics measurements for 
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air-pressure-induced surface displacements" . However, this simple model was used in the 
20th century. In order to evaluate the errors of this approach, we computed coefficients 
of linear regression between the 22 year long time series of the displacements caused by 
atmospheric pressure loading and local pressure at the station obtained by interpolation 
of the NCEP surface pressure field. Using these regression coefficients we computed a 
synthetic time series of the pressure loading for each station and estimated the admit- 
tance factors of these time series in solution G7. The results are shown in table 3. We see 
that the simple regression model works surprisingly well for vertical components. This 
test shows that the amount of the power of the signal which is present in the model, but 
is not present in the data, is increased from 5% to 12% when a regression model is used. 

Currently, there is no need to resort to a simplified linear regression. Numerical weather 
models are available on-line promptly, and a computation of a 20 years long series of 
pressure loading for all VLBI and SLR sites takes only several days at a personal computer 
using the efficient algorithm presented in the Appendix. 

Applying atmospheric pressure loading in a procedure of data reduction causes a small 
change in the resulting terrestrial reference frame: the maximum site position change 
among the stations which observed one year or longer is 2 mm, the velocity change is 
typically below 0.1 mm/yr with the maximum change of 0.4 mm/yr, and the scale factor 
is increased by 0.05 ± 0.02 ppb. Taking into account the horizontal component of the 
displacement due to atmospheric pressure loading causes rms differences in the estimates of 
polar motion and UTl at the level of 100 prad and differences in the estimates of nutation 
angles with an rms of 30 prad. The uncertainties of the EOF derived from processing 
daily VLBI sessions are currently at the level of 200-400 prad. Therefore, omission of the 



DRAFT 



February 2, 2008, 3:39am 



DRAFT 



PETROV AND BOY: ATMOSPHERIC PRESSURE LOADING X - 27 

horizontal atmospheric pressure loading is currently not a significant source of noise in 
the estimates of the EOP. 

5. Conclusions 

We found that vertical and horizontal displacements caused by atmospheric pressure 
loading currently can be computed with errors less than 15% by convolution of the surface 
pressure field from the NCEP Reanalysis model with Green's functions. Our analysis of 
VLBI observations of 40 strong stations for the time period from 1980 till 2002 demon- 
strates that on average only 5% of the power of the modeled vertical pressure loading 
signal, 16% of horizontal signal and 3% of the signal in baseline lengths is not found 
in VLBI data. Thus, all discrepancies between the observations and the model can be 
explained by known deficiencies of the model. 

Admittance factors of the time series of pressure loading vertical displacements without 
the annual constituent do not deviate from unity by more than the formal uncertainty. It 
means that at the confidence level of 5% the modeled signal is completely recovered from 
the VLBI observables. This estimate sets the upper hmit of possible errors in Green's 
functions: 4% for the radial Green's function and 12% for the horizontal Green's function. 

This allows us to conclude that on average our model of the displacements caused by 
atmospheric pressure loading is correct. At the same time for some stations the model 
does not agree with data perfectly. These stations are located either close to a coast, 
or in isolated islands, or in mountain regions. In the first case, poor modeling of the 
ocean response to atmospheric pressure forcing becomes a significant factor; in the latter 
case the spatial resolution of the weather model, 2°.5 x 2°.5, is too coarse to represent 
adequately the surface pressure field variations. 
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We have detected for the first time the horizontal component of the atmospheric pressure 
loading signal. This signal has never been before taken into account in routine reduction 
of geodetic observations. If not modeled, it adds noise to the horizontal site position with 
an rms of 0.6 mm and to the estimates of the EOP with an rms of 100 prad. According 
to the estimates of the error budget, the residual site displacements due to atmospheric 
pressure loading which is not accounted for by our model, on average have an rms of 0.4 
mm for the vertical component and 0.1 mm for the horizontal component. 

The model presented here shows a significant improvement with respect to the previous 
models. The amount of power which is present in the model, but not found in the data 
is 38% for the VDH model, 12% when the linear regression approach is used, and 5% for 
our model. 

We propose our model of the atmospheric pressure loading for use in routine processing 
of space geodesy observations. We have computed time series of the atmospheric pres- 
sure loading for all VLBI and SLR sites starting from May 1976. These time series are 
available on the Web at http://geinini.gsfc.nasa.gov/aplo and are updated daily. 
Since December 2002 the contribution of the atmospheric pressure loading displacements 
is applied on a routine basis at the VLBI analysis center of the NASA Goddard Space 
Fhght Center. 

We expect that the new numerical weather models with higher spatial and temporal 
resolution will improve the agreement of pressure loading time series with observations in 
mountainous regions. Development of ocean models forced by atmospheric pressure and 
winds will improve the pressure loading estimates for coastal and island stations. 
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Appendix A: Algorithm for the Computation of Displacements Due to 

Atmospheric pressure Loading 

We represent each component of the displacement as a sum of the contribution of the 

convolution integral over the land and over the ocean: 

u{r, t) = u^{r, t) + APo{t) Uo (Al) 
where APo{t) is the uniform sea-floor pressure and Ui^{'r,t), Uo{r) are 

n rn „ „ 

«L(r,t) = AP{f'ij,t)q{r,f'ij) cosipi j j G{ijj{f,f'ij))ds 

*=lj=l cell 

(AO.) 

n m „ „ V / 

^o(r) = q{r,f'ij) COS ifi j j G{ip{f,f"ij))ds 

and index i runs over latitude and index j runs over longitude. Here we replaced the 
integration over the sphere with a sum of integrals over small cells, q = 1 for the vertical 
component. 

Green's functions have a singularity in 0, so care must be taken in using numerical 
schemes for computing the convolution integral. Although the Green's function cannot 
be represented analytically over the whole range of its argument, we can always find a 
good approximation over a small range. We approximate the function G{iIj) ■ ip by a, 
polynomial of the third degree a + ^il^ + jil^'^ + Sip^. In order to compute the integral 
A2 over the cell, we introduce a two-dimensional Cartesian coordinate system with the 
origin in the center of the cell and the axis x towards east, the axis y towards north. We 
neglect the Earth's curvature and consider the cell as a rectangle with borders [-a, a] , [-b, 
b] on the x and y axes respectively. Then the integral of the Green's function over the 
cell with respect to a site with coordinates {xs, ys) is evaluated analytically: 
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Xi — —a — Xs 
yi = -b- Vs 



b a 



JjG{ij{xs,ys))ds= J J 

cell —b—a 



a 



+ /3 + 7vx^ + |/^ + 5{x + y ) \ dxdy 



( , 7 3\ 1 ^2 + Z22 ( , 7 1 



7 „.3^ i„ + ^21 
Xi + Zii 



+ 



, 7 3^ 1 + 2^22 / ,7 3^ 1 ?/2 + Zi2 , 

a a;2 + F ^2 Mil ^ \axi + - x{ In ^ h 

6 yi + 2;2i V 6 V yi + zu 



(2/2 - Z/i) (a;2 - 



+ I (^11 + ^22 + a^i a;2 + yi 2/2) 



^12 = v^i + y2 



Coordinates x^, ys are computed as 



2/1 



X2 \ y2 Z22 - yi Z21 J - xi[y2 zi2 - yi Zu 

X2 = a — Xs 
2/2 b-ys 



+ (A3) 



E(fy •T(f;„r) 



where T{f'ij,f) is 



|/, = iV(f',,)-f(f'i,,f) 



T{f'ij,f) 



r ,j X [r X r 
|f'y X [fx f;,-] 



(A4) 



(A5) 



and E{f'ij) , N{r'ij) are unit vectors in east and north direction with respect to the center 
of the ceU: 



Eif 




sin cos A' 
■ sin if' sin A' 
cos ip' 



(A6) 



We found that when the coefficients a;('0), /3('0), 7(V') and 5('^) are computed with the 
step 0.002 rad over the range [0, 0.16] rad, and with the step 0.02 rad over the range [0.16, 
tt], the error of the approximation of the integral A3 for a cell of size 0.044 rad (2°.5) does 
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not exceed 1%. At large angular distances we can consider the Green's function to be 
constant over the cell. For an angular distance more than 0.16 rad, taking the Green's 
function out of the integral A2 and replacing it with the value at the angular distance 
between the site and the center of the cell causes an error of less than 1%. 

Two land-sea masks are used for practical computation: coarse with the resolution of 
the surface pressure grid, and fine. If the cell of the coarse land-sea mask is completely 
land or completely sea, this cell is used for computing the integral A3. Otherwise, the 
coarse resolution cell is subdivided in smaller cells of the fine resolution grid, and the 
integral over each fine resolution cell is computed independently. The surface pressure 
is considered as defined at the corners of the coarse resolution cell. The pressure at the 
center of the cell is obtained by bi-linear interpolation. When Ui^{f) is computed, the 
cells which are over ocean are bypassed. Alternatively, the cells which are over land are 
bypassed when Uo{f) is computed. 

The computation of horizontal vectors is done separately for north and east components. 
The north and east components of the vector q{f, f') are 

qti{r.,f') = —T{f,f') ■ N{r) ^(r, r') = —T{f,f') ■ E{f) (A7) 

where T(r, r ') is defined in a way similar to A5, but with the reverse order of arguments, 
E{f),N{f) are defined according to A6, but are the unit north and east vectors for the 
site under consideration. 

Source code of the programs for computation of the displacements caused by the atmo- 
spheric pressure loading is available on the Web at http : // gemini . gsf c . nasa . gov/aplo . 
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Figure 1. a) Vertical and b) north displacement induced by atmospheric pressure loa- 
ding at the station Wettzell. Power spectrum of the c) vertical and d) north displacement. 

Figure 2. a) Vertical and b) east displacement induced by atmospheric pressure loading 
at the station Hartrao. Power spectrum of the c) vertical and d) cast displacement. 

Figure 3. Temporal autocorrelation of the vertical component of the atmospheric 
pressure loading at the station Algopark. 

Figure 4. Smoothed spatial autocorrelation of the vertical component of the atmo- 
spheric pressure loading. 

Figure 5. Map of VLBI stations used in the analysis. 

Figure 6. Histogram of the distribution of the coefficients of reduction of variances. 
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Table 1. Statistics of Displacements Caused by Atmospheric Pressure Loading and Statistics 
of the Differences in Time Series: NCEP Operational versus NCEP Reanalysis Dataset; CLIO 



Model versus IB Hypothesis. 



Station 


Lat 


Long 


Displacements 


Operationa' 


l/Reanalysis 




CLIO/IB 










r.m.s. 


rms of diff. 


correlation 


rms of diff. 


correlation 








vert. 


horz. 


vert. 


horz. 


vert. 


horz. 


vert. 


horz. 


vert. 


horz. 




Deg 


Deg 


mm 


mm 


% 


% 






% 


% 






ALGOPARK 


45°.8 


281°.9 


3.7 


0.6 


4.0 


5.7 


0.96 


0.99 


4.7 


11.0 


0.02 


0.00 


BR-VLBA 


47°.9 


240°.3 


3.2 


0.6 


8.8 


9.8 


0.96 


0.92 


6.6 


12.0 


0.02 


-0.08 


DSS45 


-35°.2 


149°.0 


2.5 


0.4 


5.8 


9.1 


0.97 


0.90 


9.0 


13.1 


0.03 


-0.13 


DSS65 


40°.2 


355°.7 


2.6 


0.6 


8.2 


12.7 


0.97 


0.77 


7.6 


15.0 


0.08 


0.02 


FD-VLBA 


30°. 5 


256°. 1 


2.4 


0.4 


10.4 


9.4 


0.93 


0.95 


9.2 


19.1 


-0.05 


-0.07 


FORTLEZA 


-3°.9 


321°.6 


1.0 


0.3 


15.5 


12.5 


0.94 


0.86 


25.4 


23.3 


0.35 


-0.07 


GILCREEK 


64°.8 


212°.5 


4.7 


0.6 


7.8 


12.3 


0.99 


0.95 


2.7 


16.2 


0.12 


0.06 


HARTRAO 


-25°.7 


27°.7 


2.3 


0.2 


7.8 


18.3 


0.88 


0.91 


10.3 


16.4 


-0.01 


0.04 


HATCREEK 


40°.6 


238°.6 


2.0 


0.5 


10.3 


9.3 


0.97 


0.94 


9.4 


14.4 


-0.04 


-0.10 


HAYSTACK 


42°.4 


288°.5 


2.8 


0.6 


4.5 


4.8 


0.92 


0.94 


11.2 


14.6 


-0.01 


0.06 


HN-VLBA 


42°. 7 


288°.0 


3.0 


0.6 


4.5 


5.2 


0.92 


0.96 


6.6 


19.6 


0.02 


0.04 


HOBART26 


-42°.6 


147°.4 


1.6 


0.4 


4.9 


8.1 


0.93 


0.95 


12.0 


33.4 


-0.05 


-0.07 


HRAS 085 


30°.5 


256°.0 


2.4 


0.4 


10.4 


9.4 


0.93 


0.95 


9.2 


19.1 


0.01 


-0.07 


KASHIMA 


35°.8 


140°.7 


1.2 


0.6 


7.4 


6.0 


0.96 


0.93 


17.3 


30.2 


0.11 


0.00 


KASHIM34 


35°.8 


140°. 7 


1.2 


0.6 


7.5 


6.0 


0.96 


0.93 


17.4 


33.8 


0.10 


0.01 


KAUAI 


22°.0 


200°.3 


0.5 


0.2 


11.7 


10.5 


0.99 


0.96 


52.4 


87.3 


-0.09 


0.06 


KOKEE 


22°.0 


200°.3 


0.5 


0.2 


11.7 


10.5 


0.99 


0.96 


52.4 


89.9 


-0.09 


0.06 


KP-VLBA 


31°.8 


248°.4 


1.9 


0.5 


13.8 


10.9 


0.87 


0.93 


12.1 


44.6 


-0.09 


-0.12 


LA-VLBA 


35°.6 


253°.7 


2.5 


0.5 


8.1 


9.6 


0.99 


0.95 


8.6 


44.1 


-0.07 


-0.06 


MATERA 


40°.5 


16°. 7 


2.3 


0.8 


9.9 


8.7 


0.92 


0.75 


9.1 


40.5 


-0.03 


-0.03 


MEDICINA 


44°3 


11°.6 


3.3 


0.9 


7.0 


9.3 


0.86 


0.75 


6.9 


38.0 


-0.01 


-0.01 


MK-VLBA 


19°.7 


204°. 5 


0.5 


0.2 


12.3 


10.6 


0.74 


0.76 


53.0 


37.5 


-0.12 


0.03 


MOJAVE12 


35°.2 


243°. 1 


1.9 


0.5 


12.4 


11.7 


0.85 


0.89 


11.4 


60.6 


-0.07 


-0.11 


NL-VLBA 


41°.6 


268°.4 


3.6 


0.7 


4.3 


6.1 


0.96 


0.94 


5.7 


44.5 


0.03 


-0.02 


NOTO 


36°.7 


15°.0 


1.4 


0.7 


12.9 


11.2 


0.98 


0.77 


13.6 


62.2 


-0.02 


-0.03 


NRAO20 


38°.2 


280°.2 


3.0 


0.5 


6.5 


7.1 


0.92 


0.97 


6.8 


58.1 


0.04 


0.03 


NRA085 3 


38°. 2 


280°.2 


3.0 


0.5 


6.5 


7.1 


0.92 


0.97 


6.8 


60.6 


0.04 


0.03 


NYALES20 


78°.9 


11°.9 


1.8 


0.7 


6.7 


7.0 


0.60 


0.83 


5.8 


62.6 


0.03 


-0.01 


ONSALA60 


57°.2 


11°.9 


4.7 


1.1 


2.8 


4.5 


0.60 


0.77 


4.7 


27.1 


-0.03 


-0.00 


OV-VLBA 


37°.0 


241°.7 


2.0 


0.5 


12.7 


10.9 


0.88 


0.91 


10.1 


28.6 


-0.06 


-0.11 


OVRO 130 


37°.0 


241°.7 


2.0 


0.5 


12.9 


10.9 


0.87 


0.91 


10.1 


33.8 


-0.06 


-0.11 


PIETOWN 


34°. 1 


251°.9 


2.3 


0.4 


9.8 


10.4 


0.90 


0.93 


8.8 


40.3 


-0.07 


-0.08 


RICHMOND 


25°.5 


279°.6 


0.9 


0.4 


7.9 


6.1 


0.81 


0.86 


26.5 


47.6 


-0.03 


0.02 


SANTIA12 


-33°.0 


289°. 3 


1.4 


0.4 


15.7 


16.0 


0.52 


0.73 


15.7 


44.6 


0.18 


-0.04 


SC-VLBA 


17°.6 


295°.4 


0.5 


0.3 


15.7 


10.0 


0.89 


0.85 


55.9 


44.6 


0.10 


0.00 


SESHAN25 


30°9 


121°.2 


3.5 


0.9 


7.6 


7.8 


0.77 


0.85 


7.1 


34.0 


0.03 


0.01 


TSUKUB32 


35°.9 


140°. 1 


1.4 


0.5 


7.0 


6.9 


0.95 


0.90 


14.4 


31.2 


0.11 


0.01 


VNDNBERG 


34°.4 


239°.4 


1.0 


0.4 


8.3 


8.8 


0.87 


0.91 


19.8 


30.2 


-0.08 


-0.12 


WESTFORD 


42°.4 


288°.5 


2.8 


0.6 


4.5 


4.8 


0.92 


0.94 


7.7 


31.7 


0.02 


0.04 


WETTZELL 


49°.0 


12°.9 


4.9 


0.9 


4.2 


6.4 


0.77 


0.80 


16.9 


28.3 


-0.03 


0.04 


Mean 






2.6 


0.6 


8.7 


9.1 


0.89 


0.89 


15.2 


39.1 


0.01 


-0.02 
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Table 2. The Global Budget of the Errors of Computation of the Atmospheric Pressure 



Loading Displacements. 



Error source 


rms 


1 J-* 

correlation 


Numerical evaluation of the integral 1% 


0.0? 


Green's function 


<2% 


± 1.0 


Land sea mask 


<5% 


0.7 


Surface pressure field 


10% 


0.9 


Ocean response 


10% 


0.0 


Total 


15% 
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Table 3. Global Admittance Factors. 



Solution 


Up 


East 


North 


Gl 


0.95 ± 0.02 


1.17 ± 0.06 


0.83 ± 0.07 


G5 


0.46 ± 0.09 


1.08 ± 0.26 


-0.89 ± 0.26 


G6 


0.98 ± 0.02 


1.20 ± 0.06 


1.02 ± 0.07 


G7 


0.88 ± 0.02 
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Table 4. V of the Estimates of Residual Amplitudes of Harmonic Site Position Variations 
without Applying the Model of Atmospheric Pressure Loading and after Applying the Model. 



Wave 


Without model 


With model 




\cy\ liorz 


vert liorz 


semi-diurnal 


2.77 1.73 


2.35 1.58 


diurnal 


4.26 2.31 


4.33 2.25 


semi-annual 


2.77 1.07 


2.31 1.10 


annual 


5.18 2.45 


6.00 2.45 
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Table 5. Admittance Factors from G2 Solution. 



Station Up East North 



ALGOPARK 


1, 


.14 


± 


0, 


.09 


1. 


,03 




0. 


18 


0. 


,82 


± 


0. 


21 


BR-VLBA 


0. 


,75 


± 


0, 


.07 


1. 


,11 


± 


0. 


13 


0. 


,42 


± 


0. 


12 


DSS45 


0, 


.85 


± 


0, 


.40 






















DSS65 


2, 


.30 


± 


0, 


.34 






















FD-VLBA 


1, 


.35 


± 


0, 


.07 


2. 


,44 


± 


0. 


,18 


1. 


,03 


± 


0. 


13 


FORTLEZA 


1 


.08 


± 


0, 


.46 






















GILCREEK 


0. 


.94 


± 


0, 


.03 


1. 


,18 


± 


0. 


11 


0. 


,35 


± 


0. 


15 


HARTRAO 


1, 


.53 


± 


0, 


.33 






















HAYSTACK 


1, 


,00 


± 


0, 


.36 






















HN-VLBA 


0. 


,12 


± 


0, 


.11 


1. 


,09 


± 


0. 


16 


1. 


,03 


± 


0. 


,18 


HRAS 085 


1 


.84 


± 


0, 


.39 






















KOKEE 


-0, 


.95 


± 


0, 


.37 






















KP-VLBA 


0. 


,92 


± 


0, 


.11 


0, 


,33 


± 


0. 


19 


0. 


,07 


± 


0. 


,14 


LA-VLBA 


0, 


.36 


± 


0, 


.06 


1. 


,33 


± 


0. 


15 


0. 


,68 


± 


0. 


12 


MATERA 


0, 


.79 


± 


0, 


.15 


0. 


,12 


± 


0. 


,28 


1. 


,17 


± 


0. 


22 


MEDICINA 


0, 


,67 


± 


0, 


.12 


0. 


,15 


± 


0. 


,25 


1. 


,20 


± 


0. 


19 


M0JAVE12 


0, 


.89 


± 


0, 


.24 


1. 


,93 


± 


0. 


39 


0. 


,23 


± 


0. 


,38 


NL-VLBA 


0, 


.91 


± 


0, 


.05 


1. 


,80 


± 


0. 


10 


0, 


,61 


± 


0. 


11 


NRAO20 


1 


.41 


± 


0, 


.08 


1. 


,99 


± 


0. 


23 


0. 


,07 


± 


0. 


21 


NRA085 3 


1 


.56 


± 


0, 


.16 


1. 


,61 


± 


0. 


42 


0. 


,67 


± 


0. 


35 


NYALES20 


0, 


.69 


± 


0, 


.12 


0. 


,60 


± 


0. 


,14 


1. 


,56 


± 


0. 


,18 


ONSALA60 


0, 


21 


± 


0, 


.08 


1. 


,03 


± 


0. 


12 


1. 


,59 


± 


0. 


13 


OV-VLBA 


-0, 


.73 


± 


0, 


.11 


0. 


,94 


± 


0. 


16 


0. 


,34 


± 


0. 


,15 


PIETOWN 


0, 


.12 


± 


0, 


.07 


1. 


,18 


± 


0. 


17 


0. 


,35 


± 


0. 


13 


SESHAN25 


4, 


.03 


± 


0, 


.39 






















TSUKUB32 


3 


.91 


± 


0, 


.38 






















WESTFORD 


1, 


.22 


± 


0, 


.06 


0. 


,68 


± 


0. 


11 


0. 


,06 


± 


0. 


13 


WETTZELL 


1. 


.14 


± 


0, 


.04 


0. 


,55 


± 


0. 


12 


0. 


,91 


± 


0. 


10 
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Fig. 1 PETROV AND BOY 
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Fig. 3 PETROV AND BOY 
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Distance (km) 
Fig. 4 PETROV AND BOY 
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